# Replication Archive for 
# Alexander Coppock, Emily Ekins and David Kirby (2018), 
# "The Long-lasting Effects of Newspaper Op-Eds on Public Opinion",
# Quarterly Journal of Political Science: Vol. 13: No. 1, pp 59-87.

rm(list = ls())

# Uncomment to install
# install.packages(c("tidyverse", "estimatr", "coefplot"))

library(tidyverse)
library(estimatr)
library(coefplot)

# Set working directory to replication archive

load("mturk_opeds_cleaned.rdata")
load("elite_opeds_cleaned.rdata")
load("dv_df.rdata")
source("opeds_source.R")

all_waves_mturk <- filter(mturk_opeds, responded_w2 & responded_w3)
all_waves_elite <- filter(elite_opeds, responded_w2)

fit_101 <- lm_robust(dv_amtrak_scale_w1 ~ as.numeric(Z == "amtrak") , data = filter(all_waves_mturk, pid_3_cat == "Democrat"))
fit_102 <- lm_robust(dv_climate_scale_w1 ~ as.numeric(Z == "climate"), data = filter(all_waves_mturk, pid_3_cat == "Democrat"))
fit_103 <- lm_robust(dv_flat_scale_w1 ~ as.numeric(Z == "paul"), data = filter(all_waves_mturk, pid_3_cat == "Democrat"))
fit_104 <- lm_robust(dv_vets_scale_w1 ~ as.numeric(Z == "veterans"), data = filter(all_waves_mturk, pid_3_cat == "Democrat"))
fit_105 <- lm_robust(dv_wall_scale_w1 ~ as.numeric(Z == "wallstreet"), data = filter(all_waves_mturk, pid_3_cat == "Democrat"))
fit_106 <- lm_robust(dv_amtrak_scale_w2 ~ as.numeric(Z == "amtrak") , data = filter(all_waves_mturk, pid_3_cat == "Democrat"))
fit_107 <- lm_robust(dv_climate_scale_w2 ~ as.numeric(Z == "climate"), data = filter(all_waves_mturk, pid_3_cat == "Democrat"))
fit_108 <- lm_robust(dv_flat_scale_w2 ~ as.numeric(Z == "paul"), data = filter(all_waves_mturk, pid_3_cat == "Democrat"))
fit_109 <- lm_robust(dv_vets_scale_w2 ~ as.numeric(Z == "veterans"), data = filter(all_waves_mturk, pid_3_cat == "Democrat"))
fit_110 <- lm_robust(dv_wall_scale_w2 ~ as.numeric(Z == "wallstreet"), data = filter(all_waves_mturk, pid_3_cat == "Democrat"))
fit_111 <- lm_robust(dv_amtrak_scale_w3 ~ as.numeric(Z == "amtrak") , data = filter(all_waves_mturk, pid_3_cat == "Democrat"))
fit_112 <- lm_robust(dv_climate_scale_w3 ~ as.numeric(Z == "climate"), data = filter(all_waves_mturk, pid_3_cat == "Democrat"))
fit_113 <- lm_robust(dv_flat_scale_w3 ~ as.numeric(Z == "paul"), data = filter(all_waves_mturk, pid_3_cat == "Democrat"))
fit_114 <- lm_robust(dv_vets_scale_w3 ~ as.numeric(Z == "veterans"), data = filter(all_waves_mturk, pid_3_cat == "Democrat"))
fit_115 <- lm_robust(dv_wall_scale_w3 ~ as.numeric(Z == "wallstreet"), data = filter(all_waves_mturk, pid_3_cat == "Democrat"))

fit_201 <- lm_robust(dv_amtrak_scale_w1 ~ as.numeric(Z == "amtrak") , data = filter(all_waves_mturk, pid_3_cat == "Republican"))
fit_202 <- lm_robust(dv_climate_scale_w1 ~ as.numeric(Z == "climate"), data = filter(all_waves_mturk, pid_3_cat == "Republican"))
fit_203 <- lm_robust(dv_flat_scale_w1 ~ as.numeric(Z == "paul"), data = filter(all_waves_mturk, pid_3_cat == "Republican"))
fit_204 <- lm_robust(dv_vets_scale_w1 ~ as.numeric(Z == "veterans"), data = filter(all_waves_mturk, pid_3_cat == "Republican"))
fit_205 <- lm_robust(dv_wall_scale_w1 ~ as.numeric(Z == "wallstreet"), data = filter(all_waves_mturk, pid_3_cat == "Republican"))
fit_206 <- lm_robust(dv_amtrak_scale_w2 ~ as.numeric(Z == "amtrak") , data = filter(all_waves_mturk, pid_3_cat == "Republican"))
fit_207 <- lm_robust(dv_climate_scale_w2 ~ as.numeric(Z == "climate"), data = filter(all_waves_mturk, pid_3_cat == "Republican"))
fit_208 <- lm_robust(dv_flat_scale_w2 ~ as.numeric(Z == "paul"), data = filter(all_waves_mturk, pid_3_cat == "Republican"))
fit_209 <- lm_robust(dv_vets_scale_w2 ~ as.numeric(Z == "veterans"), data = filter(all_waves_mturk, pid_3_cat == "Republican"))
fit_210 <- lm_robust(dv_wall_scale_w2 ~ as.numeric(Z == "wallstreet"), data = filter(all_waves_mturk, pid_3_cat == "Republican"))
fit_211 <- lm_robust(dv_amtrak_scale_w3 ~ as.numeric(Z == "amtrak") , data = filter(all_waves_mturk, pid_3_cat == "Republican"))
fit_212 <- lm_robust(dv_climate_scale_w3 ~ as.numeric(Z == "climate"), data = filter(all_waves_mturk, pid_3_cat == "Republican"))
fit_213 <- lm_robust(dv_flat_scale_w3 ~ as.numeric(Z == "paul"), data = filter(all_waves_mturk, pid_3_cat == "Republican"))
fit_214 <- lm_robust(dv_vets_scale_w3 ~ as.numeric(Z == "veterans"), data = filter(all_waves_mturk, pid_3_cat == "Republican"))
fit_215 <- lm_robust(dv_wall_scale_w3 ~ as.numeric(Z == "wallstreet"), data = filter(all_waves_mturk, pid_3_cat == "Republican"))

fit_301 <- lm_robust(dv_amtrak_scale_w1 ~ as.numeric(Z == "amtrak") , data = filter(all_waves_mturk, pid_3_cat == "Independent"))
fit_302 <- lm_robust(dv_climate_scale_w1 ~ as.numeric(Z == "climate"), data = filter(all_waves_mturk, pid_3_cat == "Independent"))
fit_303 <- lm_robust(dv_flat_scale_w1 ~ as.numeric(Z == "paul"), data = filter(all_waves_mturk, pid_3_cat == "Independent"))
fit_304 <- lm_robust(dv_vets_scale_w1 ~ as.numeric(Z == "veterans"), data = filter(all_waves_mturk, pid_3_cat == "Independent"))
fit_305 <- lm_robust(dv_wall_scale_w1 ~ as.numeric(Z == "wallstreet"), data = filter(all_waves_mturk, pid_3_cat == "Independent"))
fit_306 <- lm_robust(dv_amtrak_scale_w2 ~ as.numeric(Z == "amtrak") , data = filter(all_waves_mturk, pid_3_cat == "Independent"))
fit_307 <- lm_robust(dv_climate_scale_w2 ~ as.numeric(Z == "climate"), data = filter(all_waves_mturk, pid_3_cat == "Independent"))
fit_308 <- lm_robust(dv_flat_scale_w2 ~ as.numeric(Z == "paul"), data = filter(all_waves_mturk, pid_3_cat == "Independent"))
fit_309 <- lm_robust(dv_vets_scale_w2 ~ as.numeric(Z == "veterans"), data = filter(all_waves_mturk, pid_3_cat == "Independent"))
fit_310 <- lm_robust(dv_wall_scale_w2 ~ as.numeric(Z == "wallstreet"), data = filter(all_waves_mturk, pid_3_cat == "Independent"))
fit_311 <- lm_robust(dv_amtrak_scale_w3 ~ as.numeric(Z == "amtrak") , data = filter(all_waves_mturk, pid_3_cat == "Independent"))
fit_312 <- lm_robust(dv_climate_scale_w3 ~ as.numeric(Z == "climate"), data = filter(all_waves_mturk, pid_3_cat == "Independent"))
fit_313 <- lm_robust(dv_flat_scale_w3 ~ as.numeric(Z == "paul"), data = filter(all_waves_mturk, pid_3_cat == "Independent"))
fit_314 <- lm_robust(dv_vets_scale_w3 ~ as.numeric(Z == "veterans"), data = filter(all_waves_mturk, pid_3_cat == "Independent"))
fit_315 <- lm_robust(dv_wall_scale_w3 ~ as.numeric(Z == "wallstreet"), data = filter(all_waves_mturk, pid_3_cat == "Independent"))

fit_list_democrats_mturk <- 
  list(fit_101, fit_102, fit_103, fit_104, fit_105, 
       fit_106, fit_107, fit_108, fit_109, fit_110, 
       fit_111, fit_112, fit_113, fit_114, fit_115)
fit_list_republicans_mturk <- 
  list(fit_201, fit_202, fit_203, fit_204, fit_205, 
       fit_206, fit_207, fit_208, fit_209, fit_210, 
       fit_211, fit_212, fit_213, fit_214, fit_215)
fit_list_independents_mturk <- 
  list(fit_301, fit_302, fit_303, fit_304, fit_305, 
       fit_306, fit_307, fit_308, fit_309, fit_310, 
       fit_311, fit_312, fit_313, fit_314, fit_315)

fit_1101 <- lm_robust(dv_amtrak_scale_w1 ~ as.numeric(Z == "amtrak") , data = filter(all_waves_mturk, pid_3_cat == "Democrat"))
fit_1103 <- lm_robust(dv_flat_scale_w1 ~ as.numeric(Z == "paul"), data = filter(all_waves_mturk, pid_3_cat == "Democrat"))
fit_1104 <- lm_robust(dv_vets_scale_w1 ~ as.numeric(Z == "veterans"), data = filter(all_waves_mturk, pid_3_cat == "Democrat"))
fit_1105 <- lm_robust(dv_wall_scale_w1 ~ as.numeric(Z == "wallstreet"), data = filter(all_waves_mturk, pid_3_cat == "Democrat"))
fit_1106 <- lm_robust(dv_amtrak_scale_w2 ~ as.numeric(Z == "amtrak") , data = filter(all_waves_mturk, pid_3_cat == "Democrat"))
fit_1108 <- lm_robust(dv_flat_scale_w2 ~ as.numeric(Z == "paul"), data = filter(all_waves_mturk, pid_3_cat == "Democrat"))
fit_1109 <- lm_robust(dv_vets_scale_w2 ~ as.numeric(Z == "veterans"), data = filter(all_waves_mturk, pid_3_cat == "Democrat"))
fit_1110 <- lm_robust(dv_wall_scale_w2 ~ as.numeric(Z == "wallstreet"), data = filter(all_waves_mturk, pid_3_cat == "Democrat"))
fit_1111 <- lm_robust(dv_amtrak_scale_w3 ~ as.numeric(Z == "amtrak") , data = filter(all_waves_mturk, pid_3_cat == "Democrat"))
fit_1113 <- lm_robust(dv_flat_scale_w3 ~ as.numeric(Z == "paul"), data = filter(all_waves_mturk, pid_3_cat == "Democrat"))
fit_1114 <- lm_robust(dv_vets_scale_w3 ~ as.numeric(Z == "veterans"), data = filter(all_waves_mturk, pid_3_cat == "Democrat"))
fit_1115 <- lm_robust(dv_wall_scale_w3 ~ as.numeric(Z == "wallstreet"), data = filter(all_waves_mturk, pid_3_cat == "Democrat"))

fit_1201 <- lm_robust(dv_amtrak_scale_w1 ~ as.numeric(Z == "amtrak") , data = filter(all_waves_mturk, pid_3_cat == "Republican"))
fit_1203 <- lm_robust(dv_flat_scale_w1 ~ as.numeric(Z == "paul"), data = filter(all_waves_mturk, pid_3_cat == "Republican"))
fit_1204 <- lm_robust(dv_vets_scale_w1 ~ as.numeric(Z == "veterans"), data = filter(all_waves_mturk, pid_3_cat == "Republican"))
fit_1205 <- lm_robust(dv_wall_scale_w1 ~ as.numeric(Z == "wallstreet"), data = filter(all_waves_mturk, pid_3_cat == "Republican"))
fit_1206 <- lm_robust(dv_amtrak_scale_w2 ~ as.numeric(Z == "amtrak") , data = filter(all_waves_mturk, pid_3_cat == "Republican"))
fit_1208 <- lm_robust(dv_flat_scale_w2 ~ as.numeric(Z == "paul"), data = filter(all_waves_mturk, pid_3_cat == "Republican"))
fit_1209 <- lm_robust(dv_vets_scale_w2 ~ as.numeric(Z == "veterans"), data = filter(all_waves_mturk, pid_3_cat == "Republican"))
fit_1210 <- lm_robust(dv_wall_scale_w2 ~ as.numeric(Z == "wallstreet"), data = filter(all_waves_mturk, pid_3_cat == "Republican"))
fit_1211 <- lm_robust(dv_amtrak_scale_w3 ~ as.numeric(Z == "amtrak") , data = filter(all_waves_mturk, pid_3_cat == "Republican"))
fit_1213 <- lm_robust(dv_flat_scale_w3 ~ as.numeric(Z == "paul"), data = filter(all_waves_mturk, pid_3_cat == "Republican"))
fit_1214 <- lm_robust(dv_vets_scale_w3 ~ as.numeric(Z == "veterans"), data = filter(all_waves_mturk, pid_3_cat == "Republican"))
fit_1215 <- lm_robust(dv_wall_scale_w3 ~ as.numeric(Z == "wallstreet"), data = filter(all_waves_mturk, pid_3_cat == "Republican"))

fit_1301 <- lm_robust(dv_amtrak_scale_w1 ~ as.numeric(Z == "amtrak") , data = filter(all_waves_mturk, pid_3_cat == "Independent"))
fit_1303 <- lm_robust(dv_flat_scale_w1 ~ as.numeric(Z == "paul"), data = filter(all_waves_mturk, pid_3_cat == "Independent"))
fit_1304 <- lm_robust(dv_vets_scale_w1 ~ as.numeric(Z == "veterans"), data = filter(all_waves_mturk, pid_3_cat == "Independent"))
fit_1305 <- lm_robust(dv_wall_scale_w1 ~ as.numeric(Z == "wallstreet"), data = filter(all_waves_mturk, pid_3_cat == "Independent"))
fit_1306 <- lm_robust(dv_amtrak_scale_w2 ~ as.numeric(Z == "amtrak") , data = filter(all_waves_mturk, pid_3_cat == "Independent"))
fit_1308 <- lm_robust(dv_flat_scale_w2 ~ as.numeric(Z == "paul"), data = filter(all_waves_mturk, pid_3_cat == "Independent"))
fit_1309 <- lm_robust(dv_vets_scale_w2 ~ as.numeric(Z == "veterans"), data = filter(all_waves_mturk, pid_3_cat == "Independent"))
fit_1310 <- lm_robust(dv_wall_scale_w2 ~ as.numeric(Z == "wallstreet"), data = filter(all_waves_mturk, pid_3_cat == "Independent"))
fit_1311 <- lm_robust(dv_amtrak_scale_w3 ~ as.numeric(Z == "amtrak") , data = filter(all_waves_mturk, pid_3_cat == "Independent"))
fit_1313 <- lm_robust(dv_flat_scale_w3 ~ as.numeric(Z == "paul"), data = filter(all_waves_mturk, pid_3_cat == "Independent"))
fit_1314 <- lm_robust(dv_vets_scale_w3 ~ as.numeric(Z == "veterans"), data = filter(all_waves_mturk, pid_3_cat == "Independent"))
fit_1315 <- lm_robust(dv_wall_scale_w3 ~ as.numeric(Z == "wallstreet"), data = filter(all_waves_mturk, pid_3_cat == "Independent"))

fit_list_democrats_elite <- 
  list(fit_1101, fit_1103, fit_1104, fit_1105, 
       fit_1106, fit_1108, fit_1109, fit_1110, 
       fit_1111, fit_1113, fit_1114, fit_1115)
fit_list_republicans_elite <- 
  list(fit_1201, fit_1203, fit_1204, fit_1205, 
       fit_1206, fit_1208, fit_1209, fit_1210, 
       fit_1211, fit_1213, fit_1214, fit_1215)
fit_list_independents_elite <- 
  list(fit_1301, fit_1303, fit_1304, fit_1305, 
       fit_1306, fit_1308, fit_1309, fit_1310, 
       fit_1311, fit_1313, fit_1314, fit_1315)

by_party_df <-
  bind_rows(
    fit_list_democrats_mturk %>% map_df(tidy) %>% mutate(party = "Democrats", sample = "MTurk"),
    fit_list_republicans_mturk %>% map_df(tidy) %>% mutate(party = "Republicans", sample = "MTurk"),
    fit_list_independents_mturk %>% map_df(tidy) %>% mutate(party = "Independents", sample = "MTurk"),
    fit_list_democrats_elite %>% map_df(tidy) %>% mutate(party = "Democrats", sample = "Elite"),
    fit_list_republicans_elite %>% map_df(tidy) %>% mutate(party = "Republicans", sample = "Elite"),
    fit_list_independents_elite %>% map_df(tidy) %>% mutate(party = "Independents", sample = "Elite")
  ) %>%
  filter(term != "(Intercept)") %>%
  extract(
    outcome,
    c("junk1", "dv", "junk2", "wave"),
    "([[:alnum:]]+)_([[:alnum:]]+)_([[:alnum:]]+)_([[:alnum:]]+)"
  ) %>%
  mutate(dv = factor(
    dv,
    levels = c("amtrak", "climate", "flat", "vets", "wall"),
    labels = c("Amtrak", "Climate", "Flat Tax", "Veterans", "Wall Street")
  ),
  days = factor(wave, levels = c("w1", "w2", "w3"), labels = c(0, 10, 30)),
  days = as.numeric(as.character(days))
  )

g <- 
ggplot(by_party_df, aes(x = days, y = estimate, group = party, color = party, shape = party)) +
  scale_color_manual(values = c("blue", "purple", "red")) +
  geom_point(position = position_dodge(width = 2)) +
  geom_line() +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), 
                position = position_dodge(width = 2), width = 0) + 
  facet_grid(dv~sample, scales = "free") +
  coord_cartesian(ylim = c(0, 1)) + 
  theme_bw() +
  xlab("Days Since Treatment") +
  ylab("Average Treatment Effect Estimate") +
  theme(strip.background = element_blank(),
        legend.title = element_blank(),
        legend.position = "bottom")

g

